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Abstract 

In this article, we present a greedy algorithm based on a tensor product 
decomposition, whose aim is to compute the global minimum of a strongly 
convex energy functional. We prove the convergence of our method pro- 
vided that the gradient of the energy is Lipschitz on bounded sets. The 
main interest of this method is that it can be used for high-dimensional 
nonlinear convex problems. We illustrate this method on a prototypical 
example for uncertainty propagation on the obstacle problem. 

Keywords: Greedy algorithm; high dimension; obstacle problem; uncertainty 
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1 Introduction 

The main motivation for this work comes from two important and challenging 
problems in contemporary scientific computing: 



the uncertainty quantification for some nonlinear models in mechanics, 
and more precisely, for contact problems; 



• the computation of some high-dimensional functions in molecular dynam- 
ics. 

Concerning the first domain of application which is the main focus of this 
work, there is now a wide literature on the subject, ranging from specific ques- 
tions related to the modeling of the noise sources (in particular of their cor- 
relation), to dedicated methods for the evaluation of events with very small 
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probabilities (reliability) . The focus of this paper is rather on the development 
of methods to compute efficiently a reduced model which rapidly gives the output 
of interest as a function of the random variable which enters the input parame- 
ters, in the context of contact problems in continuum mechanics. Such a model 
can then be used to evaluate the distribution of the outputs (for a given dis- 
tribution of the input parameters), or to reduce the variance in a Monte Carlo 
computation for example. Many methods have been proposed in the literature 
to attack this problem0[¥]: stochastic collocation methods, Galerkin methods, 
perturbation methods, etc. To be more specific, let us assume that the noise 
on the parameters of the model can be modeled by a possibly large number 
of random variables T = {Ti, . . . ,Tp) G Rp, so that the quantity of interest 
(say the displacement field) u{t, x) is a function of {p + d) variables, where d is 
the dimension of the physical space. The question is then how to approximate 
a function on such a high-dimensional space. The natural idea at the basis of 
many methods is to look for the solution to this problem as a linear combination 
of tensor products: 

/ m 

w(t,x) = ^^ (7'^' 0,(0^,(2^), 

i=l j=l 

where {4'i)i<i<i and {i'j)i<j<m are bases of vector spaces of dimension I and 
m respectively which are fixed a priori, and where (J7'-')i<i<;_i<j<m are real 
numbers to be computed. This method leads to the resolution of a problem in 
a vector space of dimension N = Im which may be very large. This difficulty 
becomes all the more pregnant if p is really large, so that the solution should 
be typically approximated as a sum: 

/ / m 

u{t,x) = ^ ... ^ 5]C/^--'-^-0i^(ti)...0f^(ip)V^,(a;). (1) 

ii=l ip = lj=l 

In this case, N = IPm will be too large for a classical discretization method. 
The method we are studying is a way to circumvent this difficulty. 

The second application we have in mind is the computation of the solution 
to a high-dimensional Poisson equation arising in molecular dynamics, called 
the committor function [5]. Mathematically, this function gives the probability 
for a stochastic process to reach a given region (say ^ C M'^) before another 
one (say B C Mf^). Using Feynman-Kac formula, it can be shown that this 
function satisfies a Poisson equation in a weighted Sobolev space, with Dirichlet 
boundary conditions (namely 1 on A and on B). Typically, the stochastic 
process lives in a high-dimensional space {d is large), so that computing this 
function is a challenge. 

In both cases, the difficulty comes from the high-dimensionality of the func- 
tion to approximate. The principle of the method we are interested in is: (i) to 
rewrite the original problem as a minimization problem: 

u G argmin£(u) (2) 
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where f is a functional defined on some Hilbert space V and (ii) to expand tfie 
solution in tensor products of lower-dimensional functions 

n 

Un{t,x) ^^rk{t)sk{x). (3) 
fe=l 

In practice, for each k, the functions and Sk are computed as linear combi- 
nations of the functions of the bases ('/>i)i<,;<( and {ipj)i<j<m so that 

I 

ru{t)=Y,Km, (4) 

1=1 

and 

m 

where for each k e N*, Rk = (^fc)i<i<; e and Sk = (5'^)i<j<m G K""- 
In the end, computing the approximation ([3]) leads to a problem of dimension 
N — n{l + m) which, provided that n remains small enough, will hopefully be 
lower than the dimension of the problem obtained with the classical approach 
N — Im when the size of the bases I and m are large. 

The reduction of dimension is even more significant when we are in the case 
of equation ([I}. Indeed, the approximation ([3]) can be adapted in this case in 
the following form: 

n 

Un{t,x) = J2^kitl) ■ ■ ■ rl{tp)skix). 
fc=l 

In this case, the overall dimension of the problem will be iV = n{pl + m) instead 
of A'^ = Fm in the classical approach. 

Such a representation of a function as a sum of tensor product of other 
functions to avoid the curse of dimensionality have already been introduced in 
the literature. One approach consists in using the so-called sparse tensor product 
representation [13l [Ml [15]. If the solution u we wish to approximate is sufficiently 
regular, one does not need to use fine discretizations in each direction. This idea 
can be used for example in Galerkin-like discretizations. However, this method 
loses its efficiency in the case when the solution u is not regular enough or when 
the mesh considered is complicated. 

We adopt another approach in this article. The principle of our method is 
to determine sequentially the pairs of functions {rk, Sk) which intervene in the 
approximation ([3]) through the following minimization problem; 

{rn,Sn) G argmin 8 ?'fc(i)sfe(a;) +r{t)s{x) , (6) 

(r,s)eV,xV^ \fc=l / 



3 



where Vt and T4 in ([SJ denote respectively Hilbert spaces of functions depending 
only on the variable t or only on the variable x. 

To rewrite the two problems mentioned above as minimization problems on 
Hilbert spaces, we penalize the constraints, namely the presence of the obstacle 
for the contact problem, or the Dirichlet boundary conditions for the high- 
dimensional Poisson problem. 

The method described above has been introduced by Chinestail4 for solv- 
ing high-dimensional Fokker-Planck equations, by NouylTUj in the context of 
uncertainty quantification in mechanics, and is very much related to so-called 
greedy algorithms |1 21 [S] used in nonlinear approximation theory. The main 
contributions of this work are the following: 

• the convergence of the greedy algorithm (131)-® to the unique solution of 
^ is proved, under the key assumptions that £ is strongly convex and 
that the gradient of £ is Lispchitz on bounded sets; 

• an exponential rate of convergence is obtained in the finite dimensional 
case; 

• an adequate procedure to solve the minimization subproblem © is pro- 
posed and tested on an academic test case. 

This paper can be seen as an extension of previous works on greedy algorithms 
[T^[5] which concentrate on the linear case, namely when £{v) = — 
where || • \\v is the norm of the Hilbert space V, and L a continuous linear form 
on V. 

We would like to stress that even if all the results and proofs are provided 
in the context of tensor products of two functions, our results can be easily 
generalized to the case of tensor products of more than two functions such as 
(dl except for the results in Section 5. We have chosen not to present the results 
in this general setting for the sake of clarity. 

The paper is organized as follows. In Section 2, we introduce the general set- 
ting for the problem we consider and state the main result of this paper, namely 
the convergence of the greedy algorithm. Section 2 also presents more precisely 
the two specific examples of application we have in mind. Section 3 is devoted 
to the proof of the convergence. In Section 4, an exponential rate of convergence 
is proved, in the finite dimensional setting (i.e. when Vt and Vx are finite di- 
mensional spaces). Section 5 shows that, under specific additional assumptions 
which are typically satisfied in the context of uncertainty quantification, the 
convergence results also hold if (r„, s„) in (jG]) is only a local minimum. Finally, 
Section 6 is devoted to a discussion of the numerical implementation, as well as 
to the presentation of test cases on a toy model. 
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2 Presentation of the problem and the conver- 
gence result 

In this paper, we are interested in the convergence of a greedy algorithm for the 
minimization of higli-diniensional nonhnear convex problems. 

We first introduce the general theoretical setting in which we prove the 
convergence, then describe two prototypical examples to which our analysis can 
be applied. 

2.1 General theoretical setting 

Throughout this article, p and d denote some positive integers, and T and X 
some open sets of MP and M"^ respectively. 

Let Vt and be Hilbert spaces of real- valued functions respectively defined 
over T and X (typically or Sobolev spaces). Let and be the norms 
of Vt and V^. 

We introduce the following tensor product for all (r, s) € Vt x Vx, 



which defines a real- valued function on T x A". 
We also denote by E = {r s | {r,s) €Vt ><■ 

Let y be a Hilbert space of real-valued functions defined on T x A". The 
scalar product of V is denoted by (., .) and the associated norm by 

Let 5 be a differentiable real-valued functional defined on V. For all v & V, 
we denote by £'{v) the gradient of £ at v. 

We make the following assumptions: 

{Al) Span(S) is a dense subset of V for H-Hv^; 

{A2) for all sequences of E bounded in V, there exists a subsequence which 
weakly converges in V towards an element of S; 

(A3) the functional £ is strongly convex for ||.||y, i.e. there exists a constant 
a e IR+ for which 



The functional £ is also said to be a-convex; 

{AA) the gradient of £ is Lipschitz on bounded sets: for each bounded subset 
K of V, there exists a nonnegative constant Lk € M+ such that 




r{t)s{x) 



(7) 




(8) 



Vu,u; e V, \\£'{v)-£'{w)\\v < Lk\\v-w\\v. 



(9) 



5 



The unique global minimizer £ on V is denoted by u. Its existence and 
uniqueness are ensured by the a-convexity of the functional £: 

u = argniin£(w). 



We are going to study the following algorithm: the sequence ((r„, s„))„gpj. € 
(Vt X Fr)^ is defined recursively by 

argniin £ [y^ru ® + r ® s \ . (10) 

{r,s)eVtXV^ J 

Throughout this article, we will denote for all n e N* , 



^ru®Sk. (11) 



fc=l 

Our main result is the following theorem, whose proof is given in Section 3. 

Theorem 2.1. Under the assumptions {A\), {A2), (A3) and {AA), the itera- 
tions of the algorithm are well- defined, in the sense that 110\) has at least one 
minimizer (r„,s„). Moreover, the sequence {un)n<£N strongly converges in V 
towards u. 



Remark 2.1. For each n € N* , the minimizer of pO|) is not unique in general. In 

particular, notice that the function Vt x 14 9 (r, s) i-> f fk®Sk-\-r(^s^ 
is not convex. 



Remark 2.2. Theorem 2.1 could be generalized to the case of tensor products 
of more than two Hilbert spaces. 

Indeed, let q G N with g > 3. Let pi, • • • ,Pq be q positive integers. Let 
7i , • • ■ , 7^ be g open subsets of R^^ , • • • , R''' respectively. We consider q Hilbert 
spaces, Vi, • • • ,Vq of real- valued functions defined respectively on 71, • • • ,Tq. 
Let F be a Hilbert space of real- valued functions defined on 71 x • • • x 7^ . 
Let £■ be a real- valued differentiable functional defined on V . We denote by 
E = {r(i) ® r(«' I (r(i\ • • • ,r(«)) G x • • • x T/^} . Our algorithm can 
then easily be adapted provided that assumptions {A\), {A2), {Ai) and {A4) 

are satisfied: {r j € Vi X ■ ■ ■ X Vq are defined recursively by 

Cn-l 
^ r^k '' ■ ■ ■ -S) r'^k ^ + r'-^^ • • • (8) r^'^ 
fc=i 

Our convergence result also holds in this case. But for the sake of simplicity, 
we will limit our analysis to the case of only two Hilbert spaces. 
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Remark 2.3. Let (., .) be the scalar product defined on Span(S) as: for all 

(ri (g) 51,7-2 (g) S2) = {ri,r2)vt{si,S2)v^, 

where (., .)vt and (., denote the scalar products of Vt and Vx respectively. 
Let |j.|| be the cross-norm associated to the scalar product (.,.). The tensor 
space of Vt and Vx, denoted as ^ g) is then defined as the closure of Span(E) 
for the product norm ||.||, 

Vt^Vx = Span(I])"'". 

Let us point out that the Hilbert space V is not necessarily equal to Ft (g 14 , 
the tensor space of Vt and Vx associated to the tensor product ([7]). Indeed, an 
example where our analysis can be applied and where V Vt (E) Vx is given in 
Section 2.2.2 (see Remark 2.5.). However, the following inclusion relationship 
holds: Vt(S)Vx C V. 

Remark 2.4. If Vt and Vx are discretized in finite-dimensional spaces of dimen- 
sion I and m, our algorithm consists in solving several problems in dimension 
/ -|- m instead of solving one problem of dimension Im. Thus, we can circumvent 
the curse of high-dimensionality. 

2.2 Prototypical problems 

To prove that the general theoretical setting we described in Section 2.1 is 
satisfied on the prototypical problems we present in this section, we need the 
following lemma, which is well-known in distribution theory [11]. 

Lemma 2.1. Let U £ I?'(T x X) be a distribution such that for any functions 

([/, (p ® ^)(^v'(TxX),V(Ty:Xy) = 0- 

Then U ^ in 2?'(T x X). Moreover, for any two sequences of distribu- 
tions Rn G V {T) and Sn G V {X) such that lim„_>.oo i?„ = R in 2?'(T) and 
lim„_^oo Sn = S in V'{X), lim„^oo Rn® Sn = R® S inV {T x X). 

2.2.1 Uncertainty propagation on obstacle problems 

An example of application of our algorithm is the study of uncertainty propaga- 
tion on obstacle problems. We assume that uncertainty can be modeled by a set 
of p random variables Ti, T2, Tp, and that the random vector T = (Ti, Tp) 
takes its values in T. 

We consider also that the physical problem is defined over the domain A", 
which is supposed to be a bounded subset of M'^. If iJ is a Hilbert space of 
functions defined on A", we denote by 

L|(r, H) = {v:T^H\¥. [\\v{T)\\l] < +(X)} , 
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where E denotes the expectation with respect to the probabihty law of T, and 
\\.\\h denotes the norm of H. We endow L-^{T,H) with the scalar product 
defined by {v, w)L2^^r,H) = ^ M^), w{T))h]. 

A formulation of the obstacle problem with uncertainty is the following [7]. 
Let g G L'^{T,Hg{X)) and / £ L'^{T, H^^{X)). A membrane is stretched 
over the domain X and is deflected by some random force having pointwise 
density /(T, a:) for x d X . At the boundary dX, the membrane is fixed and in 
the interior of X the deflection is assumed to be bounded from below by the 
function g{T, x) (a random obstacle). Then the deflection z = z{T, x) is solution 
of the following obstacle problem with uncertainty (see figure 1): 

r -A^z{t,x)> f{t,x) ^ 

I z{t,x) > g{x,t) > for a.a. (t, x) e T X A", 

I iA^zit,x)+f{t,x)){z{t,x)~g{t,x))^0 J 

[ z{t,x) — ioi a.a.. {t,x) e T X dX. 

(12) 




Figure 1: Obstacle problem. 
An equivalent formulation of this problem is the following. Let us denote 

ICg = {v e L^t{T,hI{X)) I for a.a. {t,x) eTx X, v{t,x) > g{t,x)} . 
Solving the obstacle problem (jl2p consists in solving the minimization prob- 



lem 



where J'(w) 



inf J{v), 



(13) 



E 



i \VMT, x)\^dx - (/(T, .), viT, .))H-^x),Him 



One of the main difficulties of this kind of problems is their very high nonlin- 
earity. Many methods have been proposed to approximate the solution of these 
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problems in the case without uncertainty [31 [7] ■ Among them, penaHzation 
methods [71 [6] are among the most widely used. They consist in approximating 
the solution of a given obstacle problem by a sequence of solutions of penalized 
problems defined on the entire Hilbert space. 

Let p be a parameter in Such a penalized problem associated with 

problem ([T3|) may be defined as 

inf Jp{v), (14) 

v<^Ll(T.Hl[X)) 

where Jp{v) = J{v) + E [f [g{T, x) - w(T, x)]\dx\ . 

Here and below, we denote by [a]_|_ the positive part of the real number a, 
i.e. [a]+ = if a < and [a]+ = a if a > 0. 

When p goes to infinity, the solution Zp of problem ([T4|) strongly converges 
to the solution z of problem The goal of the algorithm we described in the 
previous section is to calculate the solution u — Zp oi this regularized problem 
for a given value of the parameter p. 

Let us check that the general theoretical setting we described in Section 2.1 
can be applied in this case. 

Let us consider V = L^(X,H^iX)), Vt = L^T,R), K = H^{X) and 
£iv) — Jp{v) for V & V. We have — {r ® s \ {r, s) & Vt x Vx}- We endow 
Hq{X) with the scalar product defined by (si, S2) hKx) — Jx si{x).'S/ S2{x)dx. 
In this case, we have V — Vt (E) Vx and as a consequence assumption (Al) is 

obviously satisfied. 

Besides, assumption {A2) is satisfied as well. If ((r„, s„))„gN € (Vt x Vx)^ is 
such that (||r„ (g) Sn\\v)n<£N is bounded, it is possible to extract a subsequence 
which weakly converges in V towards an element w Cz V. Besides, there exists 
a non-negative constant C £ R+ such that for all n £ N, 

\\rn(g>s„\\y = E 



X 



(r-„ «) Sn) (T, x)\^dx 
E[k„(r)p] / \WxSn{x)\^dx 



X 



< C. 



We can then choose ((r* , s* ))„gN G {Vt x VxY'^ such that (g) s* = r„ (g) s„ 
and li?'^|iL|,(r,K) = 1- The sequences (r*)„gN and (s*)„gN are then bounded 
in I/|n(T, M) and Hq{X) respectively and we can extract subsequences which 
weakly converge in L|,(T,IR) and H^{X) towards r°° € L|(r,M) and s°° € 
Hq{X) respectively. As the weak convergences in L^(T, M) and Hq{X) imply 
the convergences in the distributional sense, the sequence r* (g) s* — r„ (g) s„ 
necessarily converges towards r°°(g)s°° in 2?'(Tx A") by Lemma 2.1. As the weak 
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convergence in V also implies the convergence in the sense of the distributions, 
we obtain, by uniqueness of the hmit, w = r°° (g) s°° G S. Hence assumption 
(A2) is satisfied. 

The functional £ is differentiable and 1-convex. Indeed, for all v £ V, 
£{v) = ^\\v\\1.+ (e (/(r,.),«(T,.)>^-i(;,),^i(;,) + ^^[5(T,x)-i;(r,x)]^da; ) , 

is the sum of a 1-convex function {V 3 v ■^\\v\\y) and of a convex function 

(V^v^E [(/(T,.),^^.))^-!^,^!^ + §IA9iT,x) ~ viT,x)]ldx\). The 
functional £ therefore obeys property ([S]) with a — 1. Hence, assumption (^3) 
is satisfied. 

Let us finally check that the gradient of £ is Lipschitz. For all v,w,y G V, 



\{£\v)-£\w),y)\ < 



E 



-P 



Vj;(w(T, - w{T,x)).VxV{T,x)dx 



X 



E 



X 



{[9{T, x) - v{T, x)]+ - [g{T, x) - w{T, x)] + )y{T, x)dx 



< Wv-wWvWvWv 



-pE 



ligiT.x) - v{T,x)]+ - [giT,x) - w{T,x)] + \ \yiT,x)\ dx 



X 



For a,5 e M, it can easily be seen that |[a]+ — + | < |a — 5|. This implies 

v{T,x) - w{T,x)\ \y{T,x)\dx 



\{£'{v)-£'{w),y)\ < \\v - w\\v\\y\\v + pE 
< \\v-w\\v\\y\\v 



X 



-p{E 



\v{T,x) — w{T, x)]\'^dx 



X 



1/2 



E 



\yiT,xWdx 



X 



1/2 



The Poincare inegality in Hq{X) implies that there exists a nonnegative 
constant D E M+ such that for all h E V, 



\h{T,x)]\^dx 



X 



1/2 



< D\\h\\ 



This yields 



\{£'{v) - £'{w),y)\ < (1 + pD')\\v - w\\v\\y\\v, 



hence, 



\\£'{v) - £'{w)\\v < (1 + pD^)\\v - w\\v. 

The functional £ then obeys property ([9]) with a constant L = 1 + pD^ inde- 
pendent of the bounded set considered. 
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Thus, our obstacle problem ([T^ falls into the general theoretical setting 
introduced in Section 2.1. 

There exist several variants of the obstacle problem which could be tackled 
with our algorithm. We refer to Ref. [7J or Ref. |S1 for such examples. 

2.2.2 High-dimensional Poisson equation 

Our algorithm may also be used to calculate the solution of other problems 
than obstacle problems. Other examples are high-dimensional nonlinear Pois- 
son equations. A specific application where such high dimensional Poisson equa- 
tions arise is the calculation of the so-called committor function in molecular 
dynamics [5], which is an important quantity to compute reaction rates or to 
derive some effective dynamics for example. 

Let q G N* . The committor is the solution to the following problem: 

2; = argmini / \Wv{y)\'^ exp{-U{y)) dy 

vew ^ JR<!\(Aus) 

where q is typically large, A and B are disjoint smooth open sets of M', IJ : 
R"^ — ^ R is a given potential function such that Jj^, exp(— [/) < co and 

\ V = \ on A and u = on i? J 

For y S M'^\(^Ui?), z{y) can be interpreted as the probability that the stochastic 
process Q\ solution to 

Q\ = y-\ VC/(Q^)ds + V2W^t 
Jo 

reaches A before B. Here, Wt denotes a g-dimensional Brownian motion. 

Let p,d (z N* such that q = p + d. In this example, we consider the case 
when C = M"? \ ( A U B) is bounded. Let T and X be open convex bounded 
subsets of MP and R'^ respectively such that C C :— T x X and such that 
^{{AU B) nil) ^ where fi denotes the Lebesgue measure. We also assume 
that U e C°°(R'^). In this case, the initial problem can be rewritten as a 
minimization problem set on 

W ^ {v e H^{T X A")! w = 1 on A n r2 and w = on B n 17} , 

instead of W. Indeed, as [/ G C°°(R'^,R) and ft is bounded, there exists con- 
stants j,K > such that for all y S fi, 7 < exp{—U{y)) < k. And thus, we have 
V gW if and only if v\n € W, v\-^\^^ — 1 and wj^^fj ~ 0. 
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The penalized version of the committor problem then reads 

u = argmin £{v), (15) 

where 

" ^ ' ' ' Hy)fdy], 



£{v) = \ ( \yv{y)\^eM-U{y))dy+^ ( [ \v{y) - ifdy 
for some p > 0. 



Bnci 



Let us check that the general theoretical setting described in Section 2.1 is 
relevant for this problem. 

In this case, we consider V = H^{T x X), Vt = H^{T) and = H^{X). 
The inner products that are defined over these Hilbert spaces are the following. 
For all vi,V2 e V, ri,r2 e 14, si,S2 e V^, 



{vi,V2)v = j j {Vv\{t,x).Vv2{t,x) +vi{t,x)v2{t,x)) dtdx, 

{ri,r2)v, =^(Vri(t).Vr2(i)+ri(t)r2(i)) dt, 

{si,S2)v^ = / {'Vsi{x).Vs2{x) + Si{x)s2{x)) dx. 

Remark 2.5. Let us point out that in this case, V ^Vt®Vx. Indeed, for all 
(r, s) e 14 X 14, the V-norm of the tensor product r^s reads 

\\r (g) s\\\r = W\\'l'^{t)\\^ ^\\'l'^(x) + ll^''lli2(r)ll*lli2(Ar) + lklli2(r)ll*lli2(Ar)' 

which is not a cross-norm, equivalent to the norm induced by \\.\\vt o-nd 
over Vt^Vx, which is 

lk<S's||vt®K = Ikllvtl|s||i4- 

Indeed, let us consider T = X = (0, 1), n(i) = j sin(^^7rt) and si{x) = j sin{PiTx) 
for {t,x) e (0,1)^ and ? e N*. The sequence <H) ||y)^gf^, is bounded, but 
the sequence {\\ri\\vt\\si\\vj i^j^* is not. 

Assumption (Al) holds true, since T, = {r (g> s \ (r, s) € C°° (T) x C°° (X) } 
is such that S c S and Span (t^ is dense in H^{T x X). Hence, Span(S) is 
also dense in V. 

Let us prove that assumption {A2) also holds true. If {{rn, Sn))neN G (14 x 
VxY^ is such that (||r„ (S) s„||v')neN is bounded, we can extract a subsequence of 



12 



{rn <8) Sn)neN* which weakly converges in V towards an element w eV. Besides, 
there exists a nonnegative constant C S M+ such that for all n E N, 

\\rn<»sj^y = / (|Vr„(t)ns„(x)|2 + |r„(t)p|Vs„(a;)|2 + |r„(t)ns„(a;)|2) dtda; 
JTxX 

= ll^^n|lL2(r)ll*"lli2(A') + lkn||i2(7-)||Vs„||^2(;t.) + || »'n || 1,2 (7-) || S„ 1 1 ^2 

< c. 



We can then choose s* ))„gN G (Vf x VxY^ such that r* ®s* = r„(g)s„ and 
such that ||?'^||l2(7-) = 1. The sequences (r* )„gN and (s* )„gN are then bounded 
in iv^(T) and H^(X) and we can extract subsequences which weakly converge in 
L?{T) and H^{X) respectively towards r°° and s°° . As the weak convergences in 
L^{T) and H^{X) imply the convergences in the distributional sense, r^^® s*^ = 
Tn ® Sn necessarily converges towards r°° ® s°° in the distributional sense by 
Lemma 2.1. As the weak convergence in V also implies the convergence in the 
sense of the distributions, by uniqueness of the limit, w = r°° (g) s°°. Let us 
suppose w 0. In that case, we have r°° ^ and s°° ^ 0. Besides, we have 

Ml - ||r°°||i2(r)||V.-||i.(^) + ||r-||l,.(^)P-|li.(;,), 

hence 

IL,,I|2 

\v 



< 



II' ll/^MT) ^ ||„oo||2 

II* \\l^{X) 

As a consequence ||r°°||^i(7-) is finite and r°° e H^{T). Hence w = r°° ® s°° € 
E. If w = 0, then obviously G S. Hence, assumption {A2) holds true. 

The functional £ is differentiable and strongly convex. To prove this, it is 
sufficient to prove that there exists a constant a € such that for all v^w €V , 
{£'{v) — £'{w),v — w) > a\\v — . Indeed, there exists 7 > such that for all 
y e M', exp(— J7(y)) > 7. Thus, there exists a constant S > such that, for all 
v,w G V, 

{£'{v) - £'{w), v-w)>6( [ |V(w - w)\^ + [ \v-w\^+ [ \v - w\'- 

\Jn J Ann Jsnfi 

To prove that the functional £ is strongly convex, it is sufficient to have the 
following inequality: there exists a constant Cn £ R"^ such that for all v G 

|2 , / L.|2 ^ ^ ||„,||2 



VH'+/ kr>Cn|l«||^.(^). (16) 

n J{AuB)nn 

As T and X are bounded open convex subsets of and R*^ respectively, f2 is 
then a bounded open convex subset of R'^ such that U B) fl) 7^ and 
inequality (jl6p is a well-known Poincare-like inequality. 
Hence, assumption {A3) is satisfied. 
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Let us check that the gradient of E is Lipschitz. For all v,w,z £ V , 



\{E'{v)-E'{w),z)\ < 



< 



V(u(i, x) — w{t, x)).\/z{t, x) exp{—U{t, x)) dt dx 
{v{t, x) — w{t, x))z{t, x) dt dx 



(AuB)nn 

exp(-L/)||L^(o)||V(w - w)||i2(s:j)||Vz||i2(n) 
+p\\v - w\\L2((AuB)nn)\\z\\L2{(AuB)nn) 

< II exp(-t/)||L=e(o)||V(i; - w)||i2(f2)||Vz||i2(n) 
+p\\v - w\\L^n)\\z\\L2(^n) 

< \\v - w||y||z||y(|| exp(-[/)||i = (o) + p). 



Hence 

\\E'{v) - E'{w)\\v < (II cM-U)\\L^(n) + p)\\v - w\\v. 

The functional E therefore obeys property © with a constant i = || cxp{—U)\ \ + 
p independent of the bounded subset considered. 

Thus, the committor problem falls into the general theoretical setting intro- 
duced in Section 2.1. 



3 Proof of Theorem 2.1 

3.1 The iterations are well-defined 

We begin by proving that the iterations of the algorithm are well-defined. For 
this, we will need the following lemma. 

Lemma 3.1. Let w be a function in V . Then there exists a pair (r, s) € Vt x 
such that E{w + r <S5 s) < E(w) if and only if E'{w) ^ 0. 

Proof. Let w G V and let us suppose that f (w -I- r (g) s) > E{w) for all (r, s) € 
Vt xVx- For a given pair (r, s), for all e G M, 

E{w + er(^s)- E{w) > 0. 

As a consequence, we have the following by letting e go to 0: {E'{w), r(E}s) = 
0. This holds for all (r, s) £ Vt xV^- Hence, for all z G Span(E), we also have 
{E'{w), z) = 0, and the density of Span(E) in V , which is ensured by assumption 
{Al), yields 

E'{w) = 0. 



Conversely, let us assume that E'{w) — 0. Then, as E is a-convex, w is 
necessarily the global minimizer of E and, in particular, we have for all (r, s) G 
Vt X 14, 

E{w + r ® s)> E{w). 
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This concludes the proof. □ 

Using this lemma, the fohowing result can be proved: 

Proposition 3.1. For all n e N*, there exists a solution (r„,s„) E Vt x Vx to 
the minimization problem I110\} . 

Moreover, r„ (g) s„ ^ if and only if Un-i ^ u, where Un is defined by 

Proof. Firstly, let us prove the existence of a minimizer for problem (jlOp . 

Let new. For all (r, s) e Vt x , f (w„_i + r (g> s) > £{u). So m = 
inf eVtxv^ £{un-i +r s) exists in M. 

We then consider a minimizing sequence (r*^'^ , s^'-')/eN G {Vt x Vx f^ such that 

lim £(u„_i + r*^'^ (g) s''^) = m. 
Using ([5|) and the fact that £'(u) — 0, we have 

+ r« ® sW) - £{u) > |||m„_i + r« ® s« - 

Then the sequence (r^'^ (g) s('^);gN is bounded in F because {£{un-i + r^'^ (g) 
s*-'''))/gN is convergent and consequently bounded. 

As assumption {A2) is satisfied, we can then extract a subsequence (which 
we still denote (r*^') (g s^'^);^^) which weakly converges in V towards an element 
of S. In other words, there exist r°° S Vt and s°° € t4 such that (r*^'^ (g s'^'^);gN 
weakly converges in V towards r°° g) s°° . 

Furthermore, as the functional £ is convex and continuous on y, 

£{un-i +r°°(S) s°°) < lim £{un-i + r^'^ g) s^'^) = m. 
Hence + (g) s^) = in so that (r^, s^) is a minimizer of problem 

m- 

Let us prove now that r°° (g s°° 7^ if and only if u„_i 7^ u. 

If u„_i = u, we have f (u + r (g s) > f (u) for all (r, s) G 14 x T4 such that 
r g) s 7^ as £ is strictly convex. So a minimizer r°° (g s°° of problem (ITOl) must 
necessarily satisfy r°° ® s°° — 0. 

Conversely, if 7^ u, we have f 7^ and from Lemma 3.1, there 

exists a pair (r, s) £ Vt x Vx such that f + r (g) s) < £(w„_i). Hence, 

£{un-i + r°° g) s°°) < £(w„_i) and r°° g) s°° cannot be equal to 0. □ 

Proposition 3.2. For each n e N*, a minimizer (7'„,s„) of problem ilO\) obeys 
the following Euler equation: 

) = 0. (17) 

This result is obtained by considering the first-order conditions of the mini- 
mization problem (fTOj) . This will be useful in the proof of convergence. 
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3.2 Proof of convergence 

In this subsection, we present the different steps of the proof. 

Lemma 3.2. The series X^i^Li Ik" ^ s„||y and the sequence {£ (un)) ^^j^, are 
convergent. 

Proof. Let us set En = £{un) — £ (Sfe=i ^fc ® ^k) ■ 

Using pHl) . En < £ (u„-i + r s) for all (r, s) G Vj x 14, and in particular, 
by taking r <Si s — 0, {En)neN* is a non-increasing sequence. Moreover, it is 
bounded from below. Indeed, for all n E W , we have i?„ > £{u). Thus, it is 
convergent. 

This implies that the sequence defined as Wn — En-\ — En is nonnegative, 
converges to 0, and satisfies < +oo. 

Besides, the a-convexity of £ yields the following inequality: 

Wn > ~{£'{Un),rn «) S„) + f ||r„ (g) SnWv- 

Using the Euler equations pT)) . (£'(M„),r„ (g) s„) — 0, and thus, Wn > 
^||r„ (g) s„||^. Hence the result. □ 

Lemma 3.3. The sequence (wn)„gpj* is bounded in V . 

Proof. By a-convexity of the functional £, we have 

£{0) > £{Un)>£{u) + {£'{u),Un-u) 

1 1 1 1 2 



Thus \\u-un\\l < ^{£iO)-£{u)). 

Therefore, the sequence (M„)„gN' is bounded in V. □ 

The following estimate is essential for the proof of convergence. 

Proposition 3.3. There exists a constant A e M.^ such that, for all n E W 
and all (r, s) E Vt x Vx, 

\{£'{un^i),r(g) s)\ < A\\rn(E> Sn\\v\\r(E> s\\v. (18) 

Proof. Let M E K+ be such that for all n E N*, < M. Its existence is 

ensured by Lemma 3.3. Let N E K.+ be such that for all n E N* , ||rn(g)s„||y < N. 
Let K = B{0,M + 2N + 3) be the closed bah of V centered at of radius 
M + 2N + 3. Let L be the Lipschitz constant associated with K in 

For all (r, s) E Vt x Vx, we have £(u„_i + r g) s) - £{un-i + )>o. 

Then, by the convexity of £, we have the following inequality 

{£'{un-i + r g) s),r„ g) s„ - r g) s) < £{un-i + rn gi s„) - £{un-i + r (g s) < 0, 
which leads to 

{£'{un-i + r g) s), r g) s) > {£' { 

)• (19) 
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Let (r,s) e Vt X Vx such that ||r (g) s\\v < niax(l, |lr„ (g) s„|lv^). Wc then 
have, by using (|9]) and (fT9)) . 



-{£'{un-i),r (g) s) 



-(£'(w„_i), r (g) s) + {£'{un-i + r®s),r<^s) 
-{£'{u„-i + r <^ s),r <Si s) 



< L 

= L 

< L 
= L 

< L 

= L 



\\r (g s 
||r g) s 

||r g) s 
||r g) s 

{£'{u„ 

||r g) s 

{£'{u„ 

1 1 r g) s 



|y- (£:'(w„_i+r(gs),rg)s) 



y - {£'{un-i + r g) s), r g) s) + (£'(«„_! + r g) s), r„ g) s„) 

-1 + r g) s),r„ g) s„) 

\y - (£'(m„-i + r g) s), r„ g) s„) 

\v ~ 

) + {£'{un-i + r„ g) 

-1 + g) s„),r„ g) s„) 

y + L||r g) s - r„ g) s„||y||r„ g) s„||y 
-1 + r„ g) s„),r„ g) s„) 

y + L||r g) s - r„ g) s„||y||r„ g) s„|ly. 



The last line has been obtained by taking into account the fact that (£'(u„_i + 
rn <8) Sn),rn gi s„) = bccausc of the Euler equation ([T7]). 

Thus, for all (r, s) G x 14 such that ||r g) s||y < niax(l, ||r„ g) s„||v), 

(f'K_i),rg)s) +2.||rg)s||2, +L||r(gs||y||r„(gs„||y +L||r„(gs„||^ > 0. 

As a consequence, 

|(£'(w„_i),r g) s)| < L||rg) s||y + i||r g) s||y||r„ g) s„||y + r„ g) S„||y. 

Let (r, s) G Vf X T4 such that ||r g) s||v = 1 and t G K such that t < 
max (1, ||r„ g) s„|| y)- Then, we have 

|(f'K„i),tr(gs)| < Lt^\\r (g) sWl + Lt\\r (g) s\\ 

And, by setting t — ||r„ g) s„|| v, we obtain the following inequality for all 
(r, s) e V"f X 14 such that ||r g) s||y = 1, 

|(£'(m„_i), r(gs)\ < 3L||r„ g) s„||y||r g) s\\v. 

Of course, this inequality also holds true for all (r, s) £ 14 x Vx such that 
||r g) s||y 7^ 1. Therefore, ^ holds with A = 3L. □ 

We now state an elementary result which will be useful in the sequel. 

Lemma 3.4. Let (a„)„(=N. be a sommahle sequence ofR^. Then, there exists 
a subsequence of (na„)„gN* which converges to 0. 

Proof. If such a subsequence could not be extracted, it would imply 
3eo > , 3no G N* , Vn > no , na„ > Sq- 
Thus, the series would diverge. Hence the contradiction. □ 
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We are now in position to complete the proof of Theorem 2.1. 

Proof. By Lemma 3.2, the sequence {£{un))n£W is convergent. Let us denot its 
hmit by E. We want to prove that E = £{u). 

Firstly, for all n G W, £{un) > f (m), since u is the global minimizer of the 
functional £. By letting n go to infinity, we obtain E > £{u). 

It remains to prove that E < £{u). 

Let us first prove that {£' {un)) n^-^, weakly converges to in V. Let M e M+ 
such that for all n G N*, < Af- Its existence is ensured by Lemma 3.3. 

Let K = B{0, M + 2 + \\u\\v) be the closed bah of V centered at of radius 
M + 2+ Let L be the Lipschitz constant associated with K in Using 

([9|) and the fact that £'{u) = 0, we have ||f (it„)||y < L\\u — Un\\v and as 
(u„)„gN* is bounded in V by Lemma 3.3, we deduce that {£' {un))n<£N* is also 
bounded in V. We can then extract a subsequence of (£'(it„))„gN' which weakly 
converges in V towards w & V. By using Proposition 3.3 and by letting n go 
to infinity in (ITSl) . we deduce that {w, r (g) s) = for all (r, s) £ Vt x V^- Then, 
as Span(I]) is dense in V with assumption (Al), necessarily w = 0. Thus the 
sequence {£' {un))neN' weakly converges to in V. 

As £ is convex, we have the following inequality for all n G N* , 

£:(«„) <f(u) + (£:'(«„), w„-u)y. (20) 

Let us prove that we can extract a subsequence of ((£'(it„), M„))„gN. which 
converges to 0. Let n &N* . By using Proposition 3.3, 

n 

\{£'{un),u„)\ < ^ |(£'(u„),rfc (g) Sfc)|, 
fc=l 

n 

< ||r„+i s„+i|jy||rfc (g Sfcllv, 

k=l 

< A{n\\r„+i g) s„+i||^)i/2 \\rk g) Sk\\^ 

\k=l 

As the sequence {J2k=i Ik's ^ ^k\\v)^^j^, converges by Lemma 3.2, we have 
Sfe=i Ikfc Sk\\v ^ Sfcli Ikfc ® Sk\\v < Furthermore, we can also extract 
a subsequence from (n||r„-|_i g) Sn+i\\y)neN* which converges to (see Lemma 
3.4). 

We can then extract a subsequence from {{£' {un),Un)v)n<EN' which con- 
verges to 0. 

By letting n go to infinity in (|20p with this subsequence, we obtain that 
E < £{u). 

We have thus proved that E — £{u). 

Besides, as the functional £ is a-convex, ([8]) yields the following inequality, 

^\\u - UnWv < £{Un) - £iu), 
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which necessarily imphes that — u„||\/ converges to when n goes to infinity, 
which proves that {un)new* strongly converges towards u in V . □ 

4 Rate of convergence in the finite-dimensional 
case 

In the case when Vt and Vx are finite-dimensional, we are able to prove that the 
algorithm converges exponentially fast. 

Theorem 4.1. We assume that Vt and Vx are finite- dimensional and that as- 
sumptions (^1), (^2), (A3) and (AA) are fulfilled. Then there exist two con- 
stants T > and a G (0, 1) such that for all n G N*, 

< fK) < tct", (21) 

and 

|l,,_^„|l^<.&"/2. (22) 



Proof. Let us denote by I = dimVt and m = dim\4 • Then we can consider that 
Vt=R^,Vx= and y = K'^™ (which is implied by (Al)). 

As the spaces are finite-dimensional, all the norms are equivalent, and we 
can consider without loss of generality that ^^'^ \\-\\v are equal to 

the Frobenius norms of M', and M'^™ defined by: 

||i?||^ = R, 

\\S\\l = S^S, (23) 
Notice that for all {R, S*) e R' x M™, 

\\R<E> S\\v ^ WRS^Wlm = \\R\\l\\S\\,n- 

Let {(j)i)i<i<i and (V'jOKjXm be orthonormal bases of Vt and Vx respectively. 
Then, (0^ ® 4'j)i<i<i,i<j<m forms an orthonormal basis of V. 

Our goal is to prove that there exists a constant a G (0, 1) such that for all 

new, 

£{un) " £{u) < a {£{un-i) - £{u)) . (24) 
Let n eW . Let us notice that 

£{Un) - £{u) £{Un) - £{Un-l) + f (w„-l) - £{u). (25) 

As for all n e N*, £{un) — £{u) > 0, it is then sufficient with to prove 
that there exists A £ (0, 1) such that 

£{un) - £{un-i) < -A {£{un-i) - £{u)) , (26) 
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to have ^ with = 1 - A G (0, 1). 



Let us notice that ([8]) and (ITTt yield 

£K)-£(m„_i) < -|||r„«)s„||^. (27) 

Besides, let M £ R+ such that for all n E N* , ||u„||y < M. Its existence is 
ensured by Lemma 3.3. Let K = B{0, M + + 2) be the closed ball of V 

centered at of radius M + + 2. Let L be the Lipschitz constant of the 

gradient of £ associated to iiT in ([9|). 

Using ([9]) and the fact that £'{u) = 0, we also have, 

£{Un-l)~£iu)<L\\u-Un-l\\l. (28) 

With (|27p and ([25)) . it is sufficient to prove that there exists a constant 
K e (0, 1) such that for ah n eW, 

\\r„(S Sn\\v > K\\u-Un-l\\v, (29) 

in order to have ([26l) and hence (|24|) . 

Indeed, if ^ holds, we then have, using ^ and ([28l) . 

^(^in) - £(w„_l) < -^||r„ (g) S„||y 



As the a-convexity of £ and the fact that £'{u) = yields 

£{Un-l) - f (m) > ^11" " Un-l\\v, (30) 

inequalities ([50)1 and (|28l) then imply that ^ < L and then (|26l) holds with 
X= §;K^ e {0,1). 



Let us prove inequality (l29l) . From Proposition 3.3, estimate (|18p holds true. 
As {4>i®'4>j)i<i<i,i<j<m forms an orthonormal basis of V , we obtain, using ([T8|) . 

/ m 

ll^'("n)lly = 5]5](f'(u„),0.0V,)' 
i=l j=l 

^ m 

= ZmA^||r„+i (g) s„+i||y. 
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We then have the following estimate: 

\\S'{un}\\v < VlrnA\\rn+i «) s„+i||y. (31) 
The a-convexity of £ and estimate (|3ip lead to 



ex. 

£{Un-l)-£{u) < -{£'{Un-l),U~ Un-l) ~ -\\U~ Un-l\\v 

CX 

< VlmA\\rn <E)Sn\\v\\u- - -^H" - Un-l\\v. 



Besides, by using the fact that £(u„_i) — £{u) > 0, we obtain 

ll^-n (g) S„||y > -^=-\\u- Un-l\\v, 

2\'lmA 

which is (f29| with k = — € (0, 1) for ^ large enough. 

Hence the result. □ 

Remark 4.1. This result can be generalized to the case of tensor products of 
more than two Hilbert spaces. Indeed, with the notation of Remark 2.2, and if 
we denote li — dimVi, ■ ■ ■ ,lq — dimVq, estimate I131\) becomes 



\\£'{u.n)\\Y < ^li ■ .-IgAWrn+l «) S„+i||y, 

and the proof still holds. 

5 Case of a local minimum 

We are able to extend the results of Theorem 2.1 and Theorem 4.1 in the case 
when {rn,Sn) in (|10|) is only defined as a local minimum which ensures the 
decrease of the energy, more precisely, when (r„, s„) is defined recursively as: 

{rn, Sn) — local argminf (u„_i + r (g) s) , (32) 

such that 

£ (m„) < £ (u„_i) , (33) 

where u„ is defined as in (jlip . 

To extend these results, we will need an additional assumption (which is 
naturally fulfilled in the finite dimensional case), see Remark 5.2 below: 

{A5) There exist /3, 7 e R+ such that 

V{r,s)eVtxV., m\vM\v.<\\r^4v<j\\r\\vMv.- (34) 
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Theorem 5.1. Let us suppose that the assumptions {Al), {A2), {A3), (Ai) 
and (Ab) hold true. Then, the iterations of the algorithm described above are 
well-defined in the sense that \ 313(1 has at least one local minimizer [rn, Sn) which 
satisfies ^^) . Moreover, the sequence ( 

^")ri£N* strongly converges in V towards 

u. 

Proof. The proof is similar to the proof of Theorem 2.1 given in Section 3 except 
for Proposition 3.3 which gives estimate (IT51) : 

V(r,s) eVtX 14, \{£'{un),r(g} s)\ < yl||r„+i (g) s„+i||,/|lr ® s|ly. 

This estimate is no longer true, but we have a similar result which will be 
enough to complete the proof. Indeed, let us prove that there exists a constant 
i? e IR+ such that 

Vne N*, V(r,s) e X 14, \{£' {un),r (E> s)\ < B\\rn (E> Sn\\v\\r <E} s\\v. (35) 

Let M g M_|. such that for all n S N*, ||un||v < Its existence is ensured 
by Lemma 3.3. Let K = B{0, M + 2) be the closed ball of V centered at and 
of radius M + 2. Let L be the Lipschitz constant associated to X in ©. 

Let (r, s) e Vt X 14 and n G N* . As (r„, s„) is a local minimum of Vi x 14 3 

{y, z) ^ S (Sfe=i '''k® Sk + y® z^, there exists a constant r] £ (0, 1) such that 
for all £ G (0, 77), we have 

£ {un-i + (rn + sr) (8) (s„ + es)) > £ (u„_i + r„ (g) s„) . (36) 

Moreover, by convexity of the functional £, we have the following inequality 

£ (m„_i + (r„ + er) (g) (s„ + es)) - £ (u„-i + r„ g) s„) 

< {£'{un + s{r„ g) s + r g) s„) + £^r g) s), e(r„ g) s + r g) s„) + e^r g) s). 

(37) 

We deduce from ([55]) . ([571) and property (O that, for all e small enough so 
that ||e(r„ (g s + r (g s„) + e^r g) s||y < 1, 

< (£'(u„ + e(r„ g) s + r ® s„) + e^r g) s),e(r„ g) s + r ® s„) + £^r g) s), 

< {£'{un), e(r„ g) s + r g) s„) + e^r g) s) + L||e(r„ g) s + r g) s„) + e^r g) s||y . 

As (r„ , s„) is a local minimum of the functional 14 x I4 3 {y,z) ^ £ ^X]fe=i i^k Sk + y 'S> z 

{rn,Sn) stiU obeys the Euler equation (fTTl) and thus (£'(«„), £(r„g)s + rg)s„)) = 
0. 

Finally, we have 

e'^{£'{un),r (g) s) + ie^||r„ g) s + r (g s„ + er (g s||^ > 0. 
Dividing this expression by and letting e go to zero, we obtain 
(£'(u„), r g) s) + L||r„ (g s + r (g Sn\\v > 0, 
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which leads to 

|(f'(u„),r (g) s)| < L (||r„ ®s + r® Sn\\v + Wn ® s - r ® Sn\\v) ■ 

AU this holds without the additional assumption (p4|) for all (r, s) e V( x T4 . 
To derive estimate p5|) . we use the additional assumption we made on \\.\\v'- 

\{£'[un),r(S>s)\^''^ < VL{\\rn(E>s + r(g)Sn\\l + \\rn<E>s-r(g)Sn\\iy^^ , 

< Vl (||r„ (g)s + r(E) s„||v + ||r„ s - r (g) s„||y) , 

< 2\/I(||r„(gs||y + lir® s„||v), 

< 2/L7(||r„||yJ|s||v4 + ||r||yj|s„|lyj. 

We can then choose (r*,s,*) £ Vt x and (r*,s*) G Ft x T4 such that 
''"n ^ = '^n ® Sn and T* ® s* = T ® s and such that = ||s^||y^ < 

^^\\rn®Sn\\v and ||r*||y, = \\s*\\v^ < ^^\\r(Ss\\v Thus, 

|(rK),r®s)|i/2 = |(rK),r*®s*)|i/2, 

< 2VLj{\\r:^\vSs*\\v. + \\r*\\vAK\\vJ, 

< 4^||r„0s„||;/2||r0s|||/^ 

And in the end, we obtain estimate (|35|) with _B = 16-^j-. With this result, 
it is then possible to conclude as in the proof of Theorem 2.1. □ 

Remark 5.1. Problem (|14p falls into the scope of Theorem 5.1. On the other 
hand, this is not the case for problem (|15p . for which property (|34|) is not true 
(see Remark 2.5). We were not able to prove a similar result in the case when 
does not satisfy property ([Ml) . 

Remark 5.2. Here are two typical examples for which assumption {A5) holds 

• In the case when V = Vt g) 14, property ^S4\ l holds with /? = 7 = 1. 
This holds in uncertainty propagation problems where V = L'^{T,H) with 
H an Hilbert space of real-valued functions defined on X . Denoting by 
Vt = if.(T, M) and T4 = H, then V ^Vt®V^. 

• In other cases, to find an approximation of the global minimum of the 
energy £ , the Hilbert spaces Vt and are usually discretized in finite- 
dimensional spaces. The problem can then be rewritten as a problem over 
Vt = M} , Vx = K™ with l,m G N* , and then V is naturally defined as the 
Hilbert space V = R'^™. Then, assumptions (Al), {A2), (A3) and {A4) 
are automatically satisfied on the discretized spaces. As all the norms are 
equivalent infinite dimension, the norms onR', R™ anrfR'^™ induced by 
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the norms defined over the original Hilbert spaces Vt, Vx and V are equiva- 
lent to the Frohenius norms, defined by Ii23\) . These norms satisfy property 
since for all {R,S) e x K"\ \\RS^\\irn = Hi?!!;!!-?!!™- Hence, the 
norms induced by the norms defined on the original Hilbert spaces auto- 
matically satisfy property ^34\ l even if the property is not satisfied in the 
continuous spaces. 

As in Section 4, we can prove that the algorithm defined by ([5^ and 
converges exponentially fast in finite dimension. 

Theorem 5.2. Let us consider the algorithm defined by iS^) and i T^^j) . Let 

l,m g N*. Let Vt = M', = R" and V = K'^™. Then there exist two 
constants r > and a G (0, 1) such that for all n E N* , 

< f K) (u) < rcr", (38) 

and 

\\u-ujv<\f^<y''^^- (39) 
V a 

Proof. As the spaces are finite-dimensional, assumptions {Al), {A2), {A3), {A4) 
and {A5) are automatically fulfilled (see Remark 5.2) and estimate (I55|) holds 
true. The proof is similar to the proof of Theorem 4.1. Indeed, (l25|), ([28]) 
and still hold. Then it is sufficient to prove an inequality similar to (|29p 
to prove Theorem 5.2. However, as (1551) holds instead of (IT5|) . inequality (I5T]) is 
replaced by: 

\\£'{un)\\v <Vl^B\\r„®s„\\l, (40) 

and consequently, an inequality similar to (|29p must be obtained in another way. 

Let M S R+ such that for all n G N*, ||wn||y < Af- Its existence is ensured 
by Lemma 3.3. Let K = B{0,M -\-2 + \\u\\v) be the closed ball of V centered at 
of radius A/ + 2 + Let L be the Lipschitz constant associated with K 

in dH). 

On the one hand, using the convexity of £, © and the fact that S'{u) — 0, 
we have 

£(m„_i) - £(u„) < -(£'(u„_i),r„ (g) Sn) < L\\rn (g> s„||y||u - u„_i||y. (41) 
On the other hand, ^U^, the convexity of £, and dM]) yield 
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Then, using (gT]), we have ^ with k = , „n e (0,1). We can 

conclude as in the proof of Theorem 4.1. □ 

Remark 5.3. The results given in this section may not stand when we consider 
more than two Hilbert spaces. Indeed, the scheme of the proof of Theorem 5.1 
cannot be easily adapted and we do not necessarily have an estimate similar to 
([35]). 



6 Numerical results 

In this section, we describe how we implemented the algorithm introduced in 
Section 2 for the resolution of problem (1141) in a very simple setting, namely a 
one-dimensional membrane problem with uncertainty. We present the numerical 
results we obtained. Additional investigations to demonstrate the applicability 
and the efficiency of the procedure on high-dimensionnal problems are still re- 
quired. We however refer to Nouy[TD] for illustrations of the interest of the 
method for problems in high dimensions. 



6.1 Implementation of the algorithm 

Let us recall problem (HH). Let / e L\{T , H^^ {X)) and g e L\{r , lll{X)) . 
Let us assume that the random variable T has a probability density p(t) for 
i G T. In other words, 

P(r ev)= [ p{t) dt, 

where I? is a measurable subset of T- 

For a given value of the penalization parameter p e M+, we wish to calculate 
an approximation of the minimizer of 



where 



inf ^M, 



i ^ \V.MT. x)P dx - (/(T, .), v{T, .))H-iixiHliX) + ^Jj9iT, x) - v{T, x)]ldx 



In other words, 

|2, 



S{v) = ]-( \yMt.x)?P{t)dxdt- i {f{t,.),v{t,.))^^^^^^^^^^^^p{t)dt 
^ ' [g{t,x) — v{t,x)]'^p{t)dxdt. 



2 



XxT 



In this case, our algorithm can be rewritten in the following form. Set /o = / 
and go — g and define recursively (r„,s„) S L|,(T, M) x Hq{X) as 

irn,Sn)& argmin f„(r(g)s), 

(r,s)eL|(r,R)xHi(A') 
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with 



£n{r'»s) = ^/ \^x{r®s){t,x)\^p{t)dxdt- ( {fn^i{t,.),r(^s{t,.))H^t,x)m(x)Pii)dt 
^ Jx-kT Jr o\ ' 

+7; / [gn^i{t,x)-r'^s{t,x)]1p{t)dxdt, 
^ JxxT 

where 

fn = /n-1 + Aa;(r„ (g) S„), 

gn = gn-l-rn<^Sn- 

Indeed, 

£ (u„_i + r (g) s) = £ - ^ / [g„_i(t,x)]+p(t) dxdt + £'„(r (g) s), 

^ JxxT 

where m„ is defined as in (|lip . 

In fact, from Theorem 5.1, it is sufficient for (r„, s„) to be a local minimum 
of L|(r,M) X H^{X) 3 (r, s) 1-^ f„(r (g) s) such that: 

) < ^ / [9n-l{t,x)]lp{t)dxdt, 
^ JxxT 

which ensures ((33)l . 

We write the algorithm in the discrete case, and, for clarity, we restrict 
ourselves to the case of two open intervals T and X of R. More precisely, 
T =]t, t[ and X =]x, x[, with t, t,x,x € R, such that t<t and x<x. 

Let Z,m g N*, which will denote respectively the number of degrees of free- 
dom in the discretized spaces of Vt and Vx- Let us introduce a subdivision 
(ti)i<i<i such that t — ti < t2 < ■ ■ ■ < < ti =t and a subdivision {xj)"Li 

such that X = Xq<Xi<---< Xm < Xm+l — X. 

Let {(f)i)i<i<i C 14 and {'ipj)i<j<m C 14 be functions such that 

MU')=SH',yi<hi' <l, (42) 

and 

ipjixj,) = Sjj,, VI < j < TO, < j' < TO + 1, (43) 

and let us consider Vt — Span((/)i)i<i<; and Vx — Span(^j)i<j<m. For example, 
Pg or Qq finite elements satisfy these properties for all q gN* . 

Our goal is to find an approximation of the function u under the following 
form 

/ m 

u(t,x)«^^C/^^(/.,(i)V',(a:), (44) 
t=i j=i 
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where U = {W^)i<^<i, € M'=^'. 

Let D g M™^™ be the symmetric positive definite matrix which corresponds 
to the discretization of the one-dimensional operator —dxx in 14:, in other words, 
for all 1 < < m, 



X 



dx'>P]ix)dx^j'{x) dx. 



Let also $ = (*"')i<,,,;'<i G M'><' and * = (*^^')i<,j'<„. e 
symmetric positive definite matrices defined as 



be the 



= / Mt)^At)pit) dt, VI <i,i' < I, 



T 



and 



X 



j(x)i\} ji (x) dx, VI < j,f < m. 



With discretization the term ^ ^^-kT ^)P-P(0 '^2; is then equal 



to 



\Vxu{t,x)\^p{t) dxdt w Tr($t/i:>C/^ 



XxT 



Similarly, if we denote by F = (i^*-')i<i<i.i<j<m G M'^™ the matrix defined 
as, for all 1 < i < ^ and 1 < j < m, 

F'' = I {f{t:■):^^®M^.■))H-^(X)Ml{X)Pi^)dt: 
•J 7~ 

the term Jj- {f{t, •),u{t^ •))h-^{X) H^{A:)Pi^) approximated by 
{fit, .), u{t, .))h-^x),hI(x) Pit) dt « Tr(FC/^). 



r 



The approximation of the term ^ J^^^[g(t,x) — u(t,x)]'^p{t) dx dt is more 
subtle. Indeed, let us approximate the function g{x,t) in the discretized space 
Vt ® % by 

I m 

i=i j=i 

Given relationships and (|43l) . a natural way to define G*-' is the following 

G*-'' = g{t„ x-j) yi<i<l, l<j<m. 

It then holds 



P 



XxT 



[git, x)-u{t, x)]\p{t) dxdtK - 



XxT 



1=1 j=i 
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Let w be the function defined as w : {t,x) G TxA" i-> w{t,x) = X]jLi(^*"' ^ U^'')'t'i{t)'4'j{x) 

From assumptions (02]) and P5|) . it holds 

= [G*^' - C/*^'] + , VI < i < 1 < j < m. 

Thus, we perform a supplementary approximation, that is to approximate the 
function w itself as 

/ m 

w{t,x) sa 

/ m 



Finally, it holds 



7; [g{t,x) - u{t,x)]1p{t)dxdt « ^ 



p / I rn \ 



Tr($[G- i7] + *[G- C/]^). 



We also made the following mass lumping approximation, that is to consider 
that $ « Id; and * « Id,„. 



Then, the discretized problem ([T4|) can be rewritten as 
Find [/ e R'><™ such that 

U e argmin^,gR,x™ : - : + f [G - 1/]+ : [G - V] + , 

where for e R'^™, 

A : B = Tr(Ai?^) ^ J2 J2 ^^^^^r 

l<i<l l<j<m 

This problem is equivalent to; 

Find U e R'"""' such that UD = F + p[G - U] + . 

For each function r E Vt and s e 14, we denote by i? G R' and G R* 
their coordinates in the bases {4'i)i<i<i and {ipj)i<j<m, which are given by 

VI < i < i?^ = J r{t)Mt)dt, 
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and 

VI < j < m, — / s{x)^pj{x) dx. 
Jx 

Our algorithm can then be rewritten as: 

Choose a threshold e > and set Fq — F, Gq ^ G. At iteration n > 1: 

1. find Rn = {Rn)i<i<i and Sn — {S-li)i<j<m two vectors respectively in M} 
and M™ such that: 

(i?„,S'„)e argmin £n{R,S), 

(fl,5)GR'xR'" 

with 

£n{R,S) = i(i?5^)i?: (i?5^)-F„_i : 

+ ^[G„_i — RS^]+ : [G„_i — RS'^] + . 

2. set F„ = K_i - (i?„5,T)i?, and G„ - G„_i - i?„5^. 

3. if + /9[G,i] + || > e, proceed to iteration n + 1. Otherwise, stop. 

The remaining question is: how can we compute (i?n,5„) at step 1? This 
critical step is described in the following section. 

6.2 Computing {Rn, Sn) 
6.2.1 Fixed-point procedure 

Let us first describe a method which has been proposed by NouyflOl and Chinesta[T], 
that is the fixed-point procedure and which we use in our final numerical imple- 
mentation (see Section 6.2.2). We present this algorithm in a particular case. 
Let us consider Vt —M'', = M™ and T/ = R'^™ endowed with the Frobenius 
norms defined by We fix a given matrix M S E'^™. Let us define the 

energy functional as £{W) — \\M — W\\y for W G M'^™. In this particuler 
case, applying the greedy algorithm described above consists in computing the 
Singular Value Decomposition of the matrix M. 

In this particular case, the greedy algorithm can be rewritten in the following 
form. 

Choose a threshold e > and set Mq = M . At iteration n > 1, 

1. find two vectors i?„ and Sn respectively in M' and M™ such that 

argmin ||M„_i - . (45) 

(ii,s)effi'xM™ 

2. set M„ = M„_i - RnSl- 
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3. if ||Af„||y > e, proceed to iteration n + 1. Otherwise, stop. 
The Euler equation associated to this problem can be rewritten as 

Mn-lSn, 



1 1 Rn 
Rn 1 1 Sfi 



The method which is generally used [8^ to solve these Euler equation is a 
fixed-point algorithm, which simply reads (for a fixed n): at iteration q > 0, 
compute two vectors 5"^'+^^) e : 



such that 



||c(9)|l2 0(9+1) _ 

p(9+l)||2 0(9+1) _ NTn(9+l) 
\-tin Wy^iJn — {"'J-n-l) n.n 



(46) 



One can check [5] that this procedure is similar to the power method to com- 
pute the largest eigenvalue (and associated eigenvector) of the matrix (M„_i)-^M„_i . 

One could think of transposing this fixed-point procedure to the case of the 
obstacle problem we consider in this article. In our case, the Euler equation 

{{Rn ■ Rn)DSn = F,f_ji?„ + p[Gn-l — i?„S'J]!^-R„, 
{DSn ■ Sn)Rn = Fn-lSn + p[Gn-l — RnSn] + Sn- 

could be solved a priori with a fixed point algorithm, which, at iteration q, 
might be written as 



/r)c;'(9+i) . c(9+i)x n(9+i) _ p 0(9+1) 



:(9+i) 



(9) 



p 



Gn — 1 ~ Rn 



(9) f ^(9+1) 



j(9+l) 



Unfortunately, we have not been able to make this fully-explicit fixed point 
algorithm converge for large values of the parameter p. We therefore have 
decided to resort to a minimization procedure. 



6.2.2 Minimization procedure 



The approach we adopt then is the following. We choose an initial pair (i?„, S^) G 
and then perform a quasi-newton algorithm to find a local minimum 



pi 



of the function 
1 



(RS^) D : (RS^) - K-i : (RS^) + | [Gn-i - RS^]^ : [Gn-i - RS^]_ 
The main difficulty is to find a proper initial pair {R^n\sk°^) such that 
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< ^[^71-1]+ : [Gn-i] + , 

to ensure that the energy decreases (see ([55]) V 

Let us describe our approach in the continuous setting with the notation 
used in Section 4. It consists in finding a pair {rn \ Sn°^) € Vt x Vx such that 

f (u„_i+ri'')®sl°)) <£{u„.i). 

We notice that for (r, s) G Vt x Vx, and ij > 0, we have 

£ (u„-i + ?7r ® s) - £ (u„-i) = T] {£' (u„_i) , r (g) s) + 0(77), 

for small enough. 

The idea is then to find a pair (r, s) £ x such that (£'(u„_i), r(X)s) < 0, 
so that there exists 77 > small enough for which £ (w„_i + -qr ® s)~£ (u„_i) < 
0. Then, ri"' ® Sn^ — iqr ® s is a, good initial guess. 



Let us first consider the pair yrn \ Snj &Vt xVx such that 
si"-*) e argmin ]- \\£' (u„_i) ~r® s||^ . 



In other words, we consider {^n\sn''^ the first term of the singular value 
decomposition of £' (un-i) in V . The Euler equations then imply 



and therefore, 



By taking ri"^ ® Sn^ 
such that 







5 Sn ^ 






^g 




2 

> 0. 








V 


-Tyri"-* (g) sl^'*, there exists then 77 > small 






< 0. 



In the discrete case associated to problem (|14l) . S'i'''* j is obtained by 

taking the first term of the singular value decomposition of the matrix Fn^i + 
p[Gn-i\+- This can be done with a fixed point procedure similar to (j46p . 

Once we have this initial guess (^Rn\Sn'^^, we perform a quasi-newton 

algorithm to minimize the energy. The computations are done with the software 
Scilab[l6] and the quasi-Newton procedure is performed via the optim procedure 
of Scilab. 

Let us point out that this procedure is intrusive in general. 
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6.3 One-dimensional membrane problem 

In this section, we present the resuhs we obtained with this algorithm on the 
following membrane problem. 

We suppose A" = T = (0, 1). We consider a random variable T following a 
uniform law of probability on the interval (0, 1). We wish to study problem (|13p 
with the following values for / and 5, 

y{t,x) e (0,1)^ , f{t,x) = -1 and g{t,x) = t[sin(37ra;)]+ + (t - l)[sin(37ra;)]_. 

The negative part of a G M, i.e. [a]_ = if a > 0, and [a]_ = —a if a < 0, is 
denoted by [a]_. 

The above problem models a rope attached at x = and a; = 1 subjected 
to gravity and resting upon obstacles whose altitudes are given by g{t,x). The 
quantity u{t, x) then represents the altitude of the rope at abscissa x when 
T ^t. 

This problem is approximated by problem (jl4p with parameter p = 2500. 
The problem is discretized with a regular mesh and Pi finite elements in each 
direction. Discretization parameters are chosen as / = m = 40. 

Figure 2 represents the altitude of the obstacles given by g{t, x) for {t, x) £ 
[0,1]2. 




Figure 2: Altitude of the obstacles as a function of t and x. 



The algorithm described in the previous sections is then applied with the fol- 
lowing stopping criterion: ||F„ + p[Gn] + \\v < 10"'' with \\A\\v = ^/^r{AA'^) = 

VEtiE;=i^f, forAGM'^x'. 

Figure 3 represents the evolution of log;^o {^{un) — £{u)) and of logiQ(||F„ + 

p[Gnh\\v). 
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f terms {Rn,Sn) computed. 



We can see that our algorithm captures very quickly the main modes of the 
solution and that both the energy and the F-normof the residue ||f„+/o[G„]+||y 
converges exponentially fast, as predicted by Theorem 5.1. 

Figure 4 represents the results obtained for the solution u{t,x). Figure 5 
and 6 represents u{t, x) and g{t, x) for some special values of T. 




Figure 4: Altitude of the rope as a function of t and x. 



As we can observe, the solution does not exactly satisfies the constraint 
u{t,x) > g{t,x). This is due to the fact that we approximate a solution Up 
of the penalized problem (fT4|) for p = 2500. This is the main drawback of 
our method; we do not approximate directly the solution of the initial obstacle 
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Figure 5: Profile of u and g for T = (left) and T = 0.375 (right). 




Figure 6: Profile of u and g for T = 0.5 (left) and T = 0.625 (right). 



problem but the solution of a close regularized problem. Indeed, if we try to 
further increase the parameter p, we face the main drawback of penalization 
methods, that is the ill-conditioning of the resulting matrices. 

7 Conclusion 

In this article, we presented a greedy algorithm based on variable decomposition 
aiming at computing the global minimum of a strongly convex energy functional. 
We proved that, provided that the gradient of the energy is Lipschitz on bounded 
sets, and that the Hilbert spaces considered satisfy assumptions (Al) and {A2), 
then the approximation given by our algorithm strongly converges towards the 
desired result. One of the main advantage of the algorithm is that it can deal 
with highly nonlinear problems. We also proved that in finite dimension, this 
algorithm converges exponentially fast. 
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We applied this algorithm in the context of uncertainty quantification on 
obstacle problems. In this frame, we considered regularizations of this kind of 
problems by penalization methods. Indeed, the obstacle problem can be approx- 
imated by a global minimization problem defined on the entire Hilbert space 
of some strongly convex energy functional where the constraints of the initial 
problem are replaced by penalization terms in the expression of the functional. 
Our algorithm gives a good approximation of the solutions of the regularized 
problem. However, the problem of ill-conditioned matrices, which is inherent 
to penalization methods, limits the accuracy with which we can approach the 
solution of the initial obstacle problem. 

One way to circumvent this problem is to use augmented Lagrangian meth- 
ods (see Ref. [SI [5j [7]) instead of penalization methods. Indeed, the former 
algorithms converge towards the true solution of the initial obstacle problems. 
The adaptation of our algorithm to such methods is work in progress. 

Another extension of our work would be to consider other problems than 
obstacle problems. In Ref. [10], a similar algorithm based on Proper Generalized 
Decomposition is used to study uncertainty quantification upon a Burger type 
equation. We believe that it could be possible to extend our proof of convergence 
in the case of such hyperbolic systems. 
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